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other  four  protons  (those  from  the  four  H atoms)  form  a tetrahedral  arrangement 
around  the  N nucleus.  This  tetrahedral  arrangement  is  then  approximated  by 
a spherical  shell  at  a radius  determined  by  a variational  calculation  for 
the  total  energy  of  the  tetrahedral  ion.  (The  resulting  radius  lies  well 
out  in  the  ion  electron  cloud.)  Since  a Wigner-Seitz  polyhedra  approach  is 
used  no  specific  lattice  structure  is  considered.  We  comment  on  the  earlier 
calculation  of  Bernal  and  Massey  (BM)  and  the  more  recent  calculation  of 
Stevenson.  Emphasis  is  on  the  BM  metallic  calculation  which  we  have  essentially 
repeated.  < 

The  concept  of  exchange  is  described,  derivation  of  the  Fock  equations 
is  given,  and  the  principal  equation  of  BM  (one  electron  outside  a closed 
shell)  is  then  derived  from  the  Fock  equations. 

The  BM  equation  is  then  solved  by  computer  for  both  the  no-exchange  and 
exchange-included  cases.  Our  (unsuccessful)  attempts  to  reproduce  the  BM  re- 
sults for  the  no-exchange  case  are  described  briefly.  Our  calculation  for  the 
exchange-included  case  is  then  described  in  some  detail.  For  this  case  we  find 
the  internal  energy  at  the  equilibrium  volume  to  be  about  0.5  ev  higher  than 
that  obtained  by  BM  (our  value  is  very  close  to  Stevenson's  -5.36  ev).  Steven- 
son's equilibrium  radius  is  5.35  a.u.,  BM's  is  4.23  a.u.;  we  obtain  4.99  a.u., 
i.e.,  shifted  considerably  toward  Stevenson's  value.  In  spite  of  these  differ- 
ences in  U vs  V between  our  results  and  BM's,  our  Gibbs  energy  vs  pressure 
curve  is  about  like  theirs. 

If  the  Stevenson  mixture  curve  for  U vs  V is  accepted,  the  BM  and  the 
present  metallic  calculation  make  a transition  to  "funny  sodium"  NH/  unlikely 
(at  least  below  the  Mbar  region).  Acceptance  of  the  Stevenson  metallic  U vs  V 
curve  would  make  such  a transition  unlikely  at  any  pressure.  We  consider  tran- 
sition to  the  "funny  sodium"  form  to  be  unlikely;  transitions  to  other  metallic 
forms  are  probably  in  the  Mbar  (and  above)  range. 
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I.  INTRODUCTION 


1 

There  has  been  renewed  interest  recently  in  a possible  transition 
from  a mixture  of  ammonia  and  hydrogen  to  metallic  NH^.  Chemically, 
the  transition  can  be  expressed  by  the  reaction 

2NH3  + H2  "*  2NH4  (!) 

2 

In  a calculation  made  some  years  ago,  Bernal  and  Massey  (BM)  con- 
cluded that  such  a transition  would  occur  between  60  and  140  kbar. 
Stevenson1,  however,  concluded  that  the  transition  would  not  occur, 
at  least  not  below  1 Mbar.  (The  compact  metallic  phase  that  would 
then  form  is  unrelated  to  the  low  density  metallic  phase  considered 
here,  namely,  the  "funny  sodium"  form.) 

The  fundamentally  correct  approach  to  determine  phase  transition 
pressure  Pt  is  to  plot  Gibbs  energy  G versus  pressure  P for  each 
phase  and  define  Pt  as  the  pressure  at  which  the  two  curves  cross. 

This  follows  from  the  thermodynamic  argument  that  dG  £ 0 for  P,T 
constant  so  that  the  only  place  both  phases  can  be  in  equilibrium  is 

3 

where  dG  = 0 . A second  approach  is  to  plot  internal  energy  U vs 
volume  V for  each  phase  and  then  define  Pt  as  the  common  tangent. 

Since  dG  = 0 implies  dU  = TdS  - PdV  for  constant  P and  T,  P = -dU/dV 
if  the  TdS/dV  term  can  be  ignored,  then  the  common  slope  of  the  U 

^D.  J.  Stevenson,  Nature  258,  222  (1975). 

2M.  J.  M.  Bernal  and  H.  S.  W.  Massey,  Mon.  Not.  R.  Astr.  Soc.  114 , 

172  (1954). 

^M.  W.  Zemansky,  "Heat  and  Thermodynamics,"  McGraw-Hill  (1943),  p. 

320. 
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vs  V plot  is  the  common  pressure  of  the  G vs  P plot  and  the  two  J 

approaches  are  equivalent.  i 

A number  of  different  methods  are  frequently  utilized  to  determine  1 

i 

G (or  U) . Typically,  the  differences  between  these  various  methods  < 

consist  of  different  approximations  to  the  potentials  and/or  different 
requirements  on  the  system  wave  functions.  When  one  of  the  phases  is 
metallic  these  considerations  involve  a decision  about  the  "exchange" 

‘ 

term  (see  below),  both  whether  and  how  to  include  this  term.  Both  BM 
and  Stevenson  include  this  effect:  BM  by  using  Hartree-Fock  theory, 

i 

Stevenson  by  using  the  Kohn-Sham-Slater  "average  exchange"  which 
supposedly  includes  exchange  and  correlation  effects  in  a local  but 

self-consistent  approximation.  Since  the  latter  procedure  tends  to  1 

give  a lower  metallic  energy  than  straight  Hartree-Fock,  it  is  some- 
what curious  that  Stevenson's  NH  energy  values  are  higher  than  those 

4 

of  BM.  Both  calculations  use  a "sphericized"  NH^ + ion.  Stevenson 

4 

claims  this  causes  him  to  be  about  0.1  ev  high.  From  a rapid 

skimming  of  Chapters  10  and  11  of  ref.  5,  this  author  would  have 

★ 

expected  this  to  make  a bigger  difference.  Both  calculations  use  a 
rigid  rQ,  the  radius  of  the  smeared  out  protons  in  the  sphericized 
model.  A flexible  rQ  would  allow  lower  energy  (this  would  be  more 

^D.  J.  Stevenson  (private  communication)  (Letter  to  Curtis  Selph, 

Oct.  1976).  | 

5J.  C.  Slater,  "Quantum  Theory  of  Molecules  and  Solids,  Vol.  1,"  t 

McGraw-Hill  (1963).  i 

* Note  added  in  proof:  The  comments  in  ref.  5 are  not  applicable  here: 

Slater's  molecules  fill  the  space;  the  NH4+  ions  do  not.  Thus, 

sphericizing  makes  little  difference.  We  thank  Dr.  Stevenson  for  I 

pointing  this  out.  H 
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important  at  high  pressures  and  correcting  would  tend  to  give  G vs  P 
a less  steep  slope). ^ 

Since  Stevenson's  paper * gives  only  a bare  outline  of  his 
calculation,  whereas  the  BM  paper  gives  a relatively  complete  outline 
of  their  calculation  (including  a tabulation  of  the  effective  potential 
and  of  the  core  wave  functions)  we  have  attempted  to  repeat  the  BM 
calculation.  The  plan  of  the  present  report  is  as  follows:  In  the 
next  section  we  discuss  what  is  meant  by  "exchange".  In  section  III 
the  Fock  equations  are  derived  and  then  the  principal  equation  of  BM 
is  obtained  from  the  Fock  equations  - details  of  these  derivations 
are  given  in  the  Appendices.  Section  IV  summarizes  our  initial 
attempts  to  obtain  BM's  U vs  V curve  (without  exchange).  Section  V 
summarizes  our  work  with  exchange  included.  A summary  of  the  entire 
report  is  given  in  Section  VI. 


II.  EXCHANGE 

Exchange  is  a quantum-mechanical  term  arising  from  the  Pauli 
exclusion  principle;  there  is  no  classical  analog.  (One  can,  however, 
give  a physical  interpretation  to  the  exchange  term;  such  an  inter- 

*D.  J.  Stevenson,  Nature  258,  222  (1975). 

2 

M.  J.  M.  Bernal  and  H.  S.  W.  Massey,  Mon.  Not.  R.  Astr.  Soc. , 114, 

172  (1954). 

6Stevenson  (ref.  4)  estimates  this  correction  to  be  about  0.03  ev/ion 
i.e.,  negligible. 
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pretation  can  be  used  to  approximate  the  exchange  term  leading  to 
Slater  "local"  exchange,  etc.  Since  BM  do  not  use  this  approximation 
we  shall  not  comment  further  on  it  here.)  Since  we  are  dealing  with 
electrons  (spin  - 1/2)  the  complete  (spin  and  space  coordinates)  wave 
function  must  be  anti-symmetric.  This  is  consistent  with  the  require- 
ment that  the  probability  of  two  fermions  occupying  the  same  state 
should  be  zero. 

The  essence  of  the  origin  of  the  exchange  term  can  be  illustrated 
by  the  following  simple  example:  consider  two  electrons  with  the 
same  spin  but  represented  by  two  different  functions  f and  g.  The 
anti-symmetric  wave  function  for  the  system  is  then 

*(1,2)  = f(l)g(2)  - g(l)f (2)  (2) 

where  (1)  stands  for  the  x,y,z  coordinates  of  the  first  electron,  stc.^ 
The  coulomb  operator  part  of  /'!'*H'l'dT1dT2  gives: 

/^*(1.2)(e2/r12)'i'(l,2)dT1dT2 

= 2/g*(2)(e2/r12)g(2)dx2/f*(l)f(l)dT1  - 

2/g*(2)(e2/r12)f(2)dx2/f*(l)g(l)dT1 

where  r^2  = r^  - rt,;  dummy  variables  have  been  interchanged  as 
appropriate. 
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For  the  general,  N-particle  case,  T would  be  the  N-particle  determinant 


*(1.2. ...N) 


^(1) 

(2) 

' ' * fx(N) 

f2(l) 

f2(2) 

. , . f2(N) 

• 

9 

« 

c 

f 

fN(l) 

V2) 

. , . fN(N) 
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In  determining  ground-state  energy  via  the  variational  procedure 
we  add  Lagrange  multipliers  to  effectively  make  all  the  functions 
independent  and  vary,  say,  f*(l).  In  the  equation  for  f(l)  the  two 
terms  above  become 

2f(l)/g*(2)(e2/r12)g(2)dT2  - 2g(l)/g* (2) (e2/r12)f (2)dt2 
The  first  term  is  clearly  the  direct  coulomb  term  and  represents  the 
potential  energy  of  the  electron  associated  with  coordinate  1 in  the 
field  of  electron  2.  The  second  term  is  the  "exchange"  term.® 
Computationally,  the  significant  thing  to  notice  is  that  in  the  direct 
coulomb  term  the  entire  2-space  integral  multiplies  the  function  f 
for  which  we  are  trying  to  obtain  a solution,  whereas,  in  the  exchange 
term,  f is  inside  the  2-space  integral.  The  exchange  operator  is 
thus  a non-local  operator,  i.e.,  the  exchange  operator  in  the 
equation  for  f(l)  is  different  from  the  exchange  operator  in  the 
equation  for  g(l). 


III.  DERIVATION  OF  PRINCIPAL  EQUATIONS 
A.  Outline  of  the  Derivation  of  the  Fock  Equations 
We  desire  the  minimum  of 

E = <^|H|^/  <y|4^>  (3) 


O 

In  the  case  of  N electrons  f would  be  replaced  by  u^,  g by  Uj  and 
each  term  would  be  summed  over  all  j including  j = i.  The  exchange 
term  would  include  only  those  Uj  for  which  the  spin  is  the  same  as 
the  u^  spin.  This  last  distinction  did  not  arise  in  our  simple 
example  as  we  took  both  electrons  to  have  the  same  spin. 
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where  X2»...*xN)  = det  [4>^  (Xj)] ; Xj  contains  both  space  and  spin 

coordinates. 

H = l (-h2/2m)V?  + l U(£)  + 1/2  l Ie2/|r  - t \ (4) 

i=l  1 i=l  i=lj=l  Z 

The  three  terms  in  H are,  respectively,  kinetic,  nuclear,  and  electron- 

electron  interaction. 

It  would  thus  appear  that  the  evaluation  of  eq.  (3)  is  a form- 
idable task:  there  are  N kinetic  energy  operators,  N nuclear  operators, 
and  N(N-l)/2  electron-electron  operators;  the  determinant  of  [4>^(Xj)] 
appearing  on  the  left  and  right  of  H has  N!  terms,  each  term  containing 
N 4>(x)  factors.  Orthogonality  relations  and  recognition  of  dummy 
variables  make  the  task  tractable. 


The  orthogonality  conditions  are 

= 6km  * 

The  denominator  of  eq.  (3)  then  becomes  simply  N!.  The  first  two 
terms  in  H contain  only  one- variable  operators;  in  App.  A we  show  that 


they  give: 


and 


N 

(-*2/2m)  l ft . (rjJVj2* iCr1)dx 
j=l 

N * ^ 

^ ft j (r^UCr^^  (r1)dx1 


(A) 

(B) 


respectively;  r contains  space  coordinates  only.  The  third  term  in 
H acts  on  two  coordinates;  in  App.  B we  show  that  this  gives: 
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ie*!  J'{/  IfkC^JlMnC^ 

z k=l  £=1  •Jk  u liul2 


lrl  - r2l 


- 6 (spin  1,  spin  2)  / 


\ Crp+t  C?j)V?2>V*i> 


dTldT2^ 


rl  - r2l 


The  first  integral  in  (C)  is  called  the  "direct"  term;  the  second  is 
called  the  "exchange"  term. 

We  now  want 

6E  = 0 (6) 

6E  = dc^*  <S<*’1*  + ^2*  + ' ' ‘ + 6V  + CcomPlex  conjugate) 


These  functions  cannot  all  be  varied  independently  since  we  want 
any  variation  to  preserve  the  orthogonality  conditions  (eq.  (5)). 
Using  the  method  of  Lagrange  undetermined  multipliers  (effectively 
making  each  independent)  we  can  consider  a particular  as  an 

arbitrary  variation;  6E  = 0 then  leads  to  the  requirement  that 
(--E  -)  = 0.  This  gives 

- (fi2/2m)V124>.(?1)  + ♦ IjVjVV 

j=l  |rrr2| 


- ^6 (spin  i,  spin  j)  [/  ,^r2^i(r2^  dT2]<J>.  (^i) > = 0 (8) 

j=1  IVf2l 
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(There  is  a similar  equation  for  each  <|>^.)  We  note  that  there  is 
considerable  confusion  in  the  literature  as  to  whether  one  can 
diagonalize  the  X^j  in  the  non-closed-shell  case  (see  App.  C);  we 
shall  assume  that  we  cannot  and  proceed  to  derive  the  BM  equation 
(their  eq.  (I))  on  that  basis. 


B.  Derivation  of  Eq.  (I)  of  BM  from  the  Fock  Equations 

As  far  as  the  labelling  of  core-electrons  is  concerned  the  BM 
model  treats  metallic  NH^  as  if  it  were  "funny  sodium";  thus  the  11 
electrons  are  taken  to  be  in  a ls22s22p63s  configuration.  The  10- 
electron  ion  problem  is  first  solved  separately  (including  exchange 
terms).  The  resulting  Is,  2s,  2p  "core"  functions  are  then  taken  as 
fixed  and  an  equation  like  eq.  (8)  is  written  for  the  3s  (valence) 
electron. 

Since  we  are  writing  an  equation  for  <}>^  = ^s,  i.e.,  for  a 
function  with  £ = 0,  there  is  no  need  to  include  the  2p  term  in 
l\ij;  this  term  is  automatically  orthogonal  to  due  to  having  a 
different  £ value.  We  also  take  X = - e^.  We  substitute 

(rx)  = (Ujtr^/r^Y^te^)  (9) 

where  is  a spherical  harmonic  (a  similar  substitution  is  made  for 
4*3s (r^) ) and  write  VQ(r^)  for  Ufrj),  with 


W = - 7/rl  - 4/rc 
= - H/r, 


rl±ro 
rl  lro 


as  per  BM.  The  tilde  signifies  that  these  functions  are  not  yet  in 
atomic  units.  As  shown  in  App.  D we  obtain: 


- t 


rt2  d2u(r  ) ~ - e2z  ffCr,)  ,, 

- s -5^-  * v°(ri)u(ri)  * — r — u(ri> 


i \f,  - t 


]_— ■ — ■■ ■ YM*(e2,*2)dr2d(cos62)d*2] 


X V^W*)  - e3s“(rl)  - Xls"l5'rl’  - X2SU2s(rl> 


(10) 


It  is  legitimate  to  let  the  above  run  only  over  Is,  2s,  2p  (no  3s): 
originally  this  sum  included  ^3s  but  for  metallic  NH4  there  is  only 
one  ^3s  ^ncti°n  (just  one  electron  in  outer  shell);  thus  the  \p^s  term 
in  the  direct  coulomb  terms  is  exactly  cancelled  by  the  term  in 
the  exchange  terms. 

We  next  put  this  equation  into  atomic  units  (see  App.  E) . 


d2u(r) 


dr2 

ls,2s,2p 

-2  l 6 (spins) [ 

j 


2Z  ff(r) 

V0  (r)u(r)  + _£££ u(r) 


, UifrJufr.)  * 

/ ^.(Gp^jdrid^osepd.J.j] 

1*1  - r| 


X «j(r)Y£j(0,«  = E3su(r)  - Xlsuls(r)  - X2su2SW  • dD 

V0,  e3s,  and  the  X's  are  now  in  Rydbergs. 

The  final  step  is  to  manipulate  the  exchange  term  to  obtain  a 
purely  radial  equation.  This  is  done  in  App.  F.  A key  point  here  is 
that  we  are  dealing  with  closed  H subshells;  this  makes  (^i *^1) X 

Y^j (0 , 4>)  invariant  under  rotations,  i.e.,  dependent  only  on  the  angle 
between  r and  r^.  We  obtain: 
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9V.  Fock,  Z.  f.  Phys.  81^,  195  (1933). 

°A  number  of  other  references  have  eq.  (12)  as  we've  written  it 
(without  the  (2£  + 1)  factor  in  the  exchange  term).  This  would 
appear  to  be  a typographic  error  in  the  BM  paper:  For  an  a value  of 
3.68  a.u.,  BM  get  a ground  state  E (including  exchange)  » -7.5  ev. 

We  obtain  -6.90  ev.  using  the  correct  expression  and  -9.65  ev.  using  » 

the  (2JL  +1)  factor.  The  much  smaller  discrepancy  between  our 
(correct  factor)  value  and  their  value  indicates  that  they  probably 
used  the  correct  factor  in  their  actual  computation. 
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IV.  RESULTS  - OUR  CALCULATION  IGNORING  EXCHANGE 


In  this  section  we  describe  our  (unsuccessful)  attempts  to 
reproduce  BM's  curve  I:  ground-state  energy  vs  "a"  (the  Wigner- 
Seitz  radius)  with  exchange  ignored. 

A.  Calculation  with  the  Lagrange  Multipliers  A^s  and  X2s  Set  Equal 
to  Zero 

Without  exchange  and  with  Als  and  X2s  set  (arbitrarily)  to  zero. 


eq.  (I)  of  BM  is: 


[-d2/dr2  + V(r)]u(r)  = eu(r)  , 

(13) 

with  boundary  conditions: 

u(0)  - 0 

(14a) 

du(r)/dr|a  = 0 

(14b) 

The  computer  program  "WAVEB"  is  used  to  solve  eq.  (13).  This 
program,  based  on  the  Numerov  technique,  is  outlined  in  App.  G. 
Results**  for  the  two  values  of  "a"  used  are  shown  in  Table  1. 


TABLE  1.  GROUND-STATE  ENERGIES  WITHOUT  EXCHANGE 
Xls  = X2s  = 0 

a e(ev.)  e(ev.) 

(a.u.) This  calc. BM 

3.68  2.74  -1.26 

4.80  -2.11  -4.37 
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All  BM  energy  values  listed  in  this  report  are  estimated  from  their 
Figure  1. 
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/ 


Our  Ac  between  the  two  "a"  values  is  4.85  ev,  not  even  close  to  the  BM 
value  of  3.11  ev.  (see  their  Figure  1,  curve  I). 

B.  "Expectation"  Value 

Our  next  step  was  to  orthogonal ize  the  u(r)  determined  in 

section  A to  uls(r)  and  u2s(r)  using  a Schmidt  procedure  and  to  then 

find  the  "expectation"  value: 

e = ^■‘■(r)  | -d2/dr2  + V(r)|ui(r)^>  (15) 

12 

The  computer  program  "DSPLSW",  based  on  a spline  routine  is  used 

13 

to  evaluate  the  numerical  derivative  involved  here  . The  computer 
program  "EINL"  evaluates  eq.  (15).  Three  comments  on  these  programs 
are  in  order: 

(1)  An  attempt  to  avoid  taking  any  derivative  through  expansion 
of  u -1- (r)  in  terms  of  u,  Ujs,  and  u2s  would  involve  considerable 
algebra  plus  values  for  elg  and  e2s* 

(2)  In  using  DSPLSW,  dux  (r)/dr|0  is  found  graphically. 

(3)  The  u2s(r),  written  as  P2g(r)  2n  notati°n>  is  not 
orthogonal  to  their  ulg(r).  It  is  straightforward  but  tedious  to 
show  that  one  obtains  the  same  u (r)  whether  or  not  one  first 
orthogonalizes  u2g  to  u^s. 

i 2 

We  are  indebted  to  Ray  Scanlon  for  furnishing  this  program  and  for 
discussions  related  to  its  use. 
i ^Integration  by  parts  allows  us  to  write 

/ uJ-(d2u4’/dr2)dr  = u^du-Vdr)  | - / (du-Vdr)2dr  . 

o 0 o 

The  boundary  conditions  u-KO)  = 0,  (duA/dr)|a  = 0,  make  the  first  term 
vanish;  thus,  only  the  first  derivative  of  uA  is  needed. 
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Results  for  the  same  two  "a"  values  employed  in  Section  A,  above, 
are  shown  in  Table  2. 

TABLE  2.  <(e]>  WITHOUT  EXCHANGE 

(See  Eq.  (15)  of  the  Text) 

a <e>  (ev.) 

(a.u.) This  Calc. 

3.68  1.11 

4.80  -2.71 

Our  Ae  between  the  two  "a"  values  is  approximately  3.8  ev,  again  not  very 
close  to  BM's  Ae  of  3.1  ev. 

C.  No-Exchange  Calculation  Including  Xls  and  X2s 

From  the  failure  of  the  two  previous  approaches  to  reproduce  the 
BM  no-exchange-term  curve,  it  appeared  one  should  include  X^s  and  X2s 
and  use  a "3-dimensional"  search  technique  (for  X^,  and  e) ; 

i.e.,  we  are  now  solving: 

[-d2/dr2  + V (r) ]u(r)  = eu(r)  - *lsuls(r)  - X2su2s(r)  (16) 
with  conditions 

(ululs)  = <Hu2}>  = 0 (17) 

in  addition  to  the  boundary  conditions  expressed  by  eqs.  (14a)  and 

(14b).  We  use  the  computer  program  "TWAVE".  The  bulk  of  the  3- 
dimensional  search  logic  is  contained  in  a subroutine  utilizing  a 
"Rosenbrock"  technique.1**  A narrative  description  of  the  TWAVE 
program  is  given  in  App.  H. 

1^The  Rosenbrock-search  subroutine  used  here  was  written  by  Ray  Scanlon. 
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Our  Ae  is  4.25  ev.,  not  very  near  BM's  Ae  of  3.11  ev.  At  this  time 
we  have  no  explanation  as  to  why  our  Ae  as  calculated  here  differs  so 
much  from  that  of  BM. 

For  completeness,  ant  since  the  u(r)  of  eq.  (16)  is  needed  as 
starting  input  for  the  iterative  solution  including  exchange 
(discussed  in  the  next  section)  we've  calculated  e without  exchange 
for  a number  of  values  of  "a."  These  results  and  those  of  BM  are 
listed  in  columns  three  and  four,  respectively,  of  Table  4;  our  no- 
exchange results  are  also  compared  with  those  of  BM  in  the  upper  part 
of  Figure  1. 


TABLE  4.  AMMONIUM  RESULTS 


a 

a.u. ) 

Vol/N  atom 
(10-^^cm^) 

No  exchange 
This  calc. 

(ev.) 

BM 

With  Exchange 
This  calc. 

(ev.) 

BM 

2.24 

6.98 

38.82 

>10. 

2.610 

>3. 

2.56 

10.41 

22.23 

>7. 

-3.226 

~1.0 

2.88 

14.83 

12.60 

>5. 

-5.573 

-4.20 

3.36 

23.55 

4.670 

0.95 

-6.721 

-7.05 

3.68 

30.93 

1.619 

-1.26 

-6.898 

-7.52 

4.00 

39.72 

- .337 

-2.63 

-6.886 

-7.60 

4.48 

55.81 

-1.830 

-3.95 

-6.717 

-7.32 

4.80 

68.64 

-2.627  | 

-4.37 

-6.584 

-6.95 

5.12 

83.31 

-3.136 

-4.70 

-6.431 

-6.65 

5.44 

100.0 

-3.468 

-4.75 

-6.282 

-6.35 

V.  OUR  CALCULATION  INCLUDING  EXCHANGE 
A.  Procedure 

In  this  section  we  are  concerned  with  the  solution  of  eq.  (12) 
with  the  right-hand  side  included;  u(r)  must  also  satisfy  the  boundary 
conditions  expressed  by  eqs.  (14a)  and  (14b).  We  follow  the  procedure 
as  outlined  by  BM.  To  start  the  iterative  process  we  first  use  the 
solution  obtained  ignoring  exchange  for  Uj-0j(r),  namely,  the  solution 
to  eq.  (16)  obtained  using  TWAVE.  The  u(r)  obtained  using  TWAVE  is 
first  orthogonalized  to  Ujs(r)  and  U2g(r)  and  then  normalized.  This 
uf(o) (r)  function  is  then  used  to  calculate  the  exchange  term  g(r). 

We  then  solve  eq.  (12)  with  Xls  = X2S  = 0* 
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(18) 


with 


[-d2/ dr2  + V(r)  - e]ur,  ..(r)  = -g,._  (r) 


|U(k+l) 


’(k) 


Is , 2s , 2p  r n 

g(k)(r)  = “ I 2 {r"  'V 


nJl 


JG  ' 00 


+ rVaunJt(x)uJ-(;k)(x)x"i"ldx}unil(r)  (19) 

The  u(r)  found  by  solving  eq.  (18)  is  then  orthogonal ized  to  ul5(r) 
and  U2S(r)  and  then  renormalized.  This  u-1-^  (r)  is  substituted  into 
eq.  (19)  and  the  resulting  g^fr)  is  used  as  input  to  eq.  (18)  to 
obtain  u^+i)(r)-  This  process  is  continued  until  "convergence"  is 
obtained.  We  have  used  four  cycles  in  most  cases;  convergence 
appears  to  be  obtained  in  three  cycles.  Convergence  here  is  in  the 
sense  that  the  orthogonal  ized  and  normalized  u-l-^+i)  (r)  s u-'-^fr) 
and  the  unorthogonalized  u^k+1-j(r)  a u^fr).  The  (unorthogonalized) 
u(r)  solution  of  eq.  (18)  does  not  appear  to  be  converging  to  ul- (r) , 
i.e.,  orthogonalization  remains  significant  no  matter  how  large  k 
becomes. 

The  operations  above  are  performed  in  the  computer  program 
"COMBINE":  orthogonalization  and  renormalization  are  done  in  sub- 
routine RENOM  using  a Schmidt  procedure;  g(r)  as  given  in  Eq.  (19) 
is  evaluated  in  subroutine  XCHNG;  eq.  (18)  is  solved  in  MAIN  of 
COMBINE  using  a Numerov  technique  similar  to  that  of  WAVEB  (see 
App.  G). 


15It  is  not  at  all  clear  that  this  procedure  of  solving  eq.  (12)  with 
Xis  = \2s  = 0 and  then  orthogonal  zing  the  resultant  u(r)  to  u^s  and 
U2S  is  valid;  we  comment  on  this  at  some  length  in  App.  I. 
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B.  Results  (Including  Exchange) 


Columns  five  and  six  of  Table  4 give  our  results  and  those  of  BM 
for  various  values  of  "a"  for  the  exchanged-included  case;  these  results 
are  also  compared  in  the  lower  part  of  Figure  1.  Our  values  are 
roughly  0.5  ev.  above  BM's  over  most  of  the  important  part  of  the 
volume  range  (20  to  60  on  the  horizontal  scale  of  Figure  1). 

BM  calculated  total  energy  by  adding  the  mean  Fermi  energy  Gp 
to  the  ground-state  energy  including  exchange.  The  mean  Fermi  energy 
was  calculated  on  a free-electron  basis  following  Mott.^  Since  our 
ground-state  curve  is  some  0.5  ev.  higher  than  BM's  and  since  our  Gp 
is  identical  to  theirs,  our  total  (internal)  energy  U will  also  be 
some  0.5  ev.  above  theirs.  Our  values  for  U are  compared  with  those 
of  BM  in  Table  5. 

To  obtain  the  difference  in  energy  between  metallic  NH4  and  the 
1 

NHj  ~ ~2  ^2  mixture  we  follow  BM  (see  their  eq.  (3)). 

[NH4,  metal]  = [NH^,  gas]  - Ap  + (g  + Gp)0  (20a) 

where  Ap  is  the  proton  affinity  of  NHj. 

1 1 
[NH3,  crystal]  + 2 [H2,  crystal]  = [NHj,  gas]  - IH  - 2 D^H2^  “ B 

(20b) 

where  1^  is  the  ionization  energy  of  atomic  hydrogen,  D(H2)  is  the 
dissociation  energy  of  H2  and 

1 

B = B . E . (NHj , crystal)  + J B.E.(H2,  crystal)  , (20c) 

l^N.F.  Mott  and  H.  Jones,  "Theory  of  the  Properties  of  Metals  and 
Alloys,"  Oxford  (1936).  (Dover  (1958),  pp.  54,55). 
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TABLE  5.  AMMONIUM  RESULTS 
U vs  V 


Vol/N  atom 
(10-24  cm3) 

eF 

(ev) 

Total  energy 
This  calc 

, U (ev) 

BM 

6.98 

5.963 

8.573 

>4. 

10.41 

4.565 

1.339 

>2. 

14.83 

3.607 

-1.966 

-0.59 

23.55 

2 650 

-4.071 

-4.40 

30.93 

2.209 

-4.689 

-5.31 

39.72 

1.870 

-5.016 

-5.73 

55.81 

1.491 

-5.226 

-5.83 

68.64 

1.299 

-5.285 

-5.65 

83.31 

1.141 

-5.290 

-5.51 

100.0 

1.011 

-5.271 

-5.45 

(3/5)Emax 

(3/5)  (3/Tr)2/3(TT2'h2/2m)  (N/Q)2/3, 

using  eq.  (19), 

Chap  II,  : 

N/ft  = (#  of  electrons)/ (unit  volume)  = 3/4ira3 
2/3 

Ep  = (3/10)  (9tr/4)  (ft2/ma2),  the  expression  used  by  BM. 

Then,  ef  = (1.10) (27.2)/a2  = (29.92)/a2 
U = Ground  state  energy  (including  exchange)  + ep 


16N.  F.  Mott  and  H.  Jones, 
Alloys,"  Oxford  (1936). 


"Theory  of  the  Properties  of  Metals  and 
(Dover  (1958),  pp.  54,55). 
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the  sum  of  the  respective  binding  energies.  The  desired  mixture 
energy  is  then 

1 

[mixture]  = [NHj,  gas]  - IH  - J D(H2)  - B - e12  (20d) 

where  the  positive  quantity  e ^ is  the  reduction  in  energy  (per  NH3 

molecule)  due  to  mixing  the  two  molecular  crystals  (our  notation 

differs  somewhat  from  BM's).  The  various  entities  in  eqs.  (20a)  - 

(20d)  are  depicted  in  Figure  2;  this  picture  is  consistent  with  eq. 

(3)  of  Stevenson*  (Stevenson's  B includes  our  e^)- 

Following  BM,  we  take  the  various  binding  energies  (B  and  e^) 

as  being  negligible.  Eqs.  (20a)  and  (20d)  then  give: 

[NH4,  metal]  - [mixture]  = (e  + Ep)0  - Ap  + + °/2  (21) 

BM  obtained  a minimum  U for  the  metal  some  0.86  ev.  above  the  minimum 

17 

for  the  mixture;  we  obtain  ~ 1.4  ev.  for  this  difference.  Our  total 

(internal)  energy  curve  is  compared  with  BM's  in  Figure  3.  (Stevenson's 

metallic  and  mixture  curves  are  also  shown.)  The  zero  of  energy  for 

Figure  3 is  that  of  infinitely  dispersed  NH4+  ions  and  electrons: 

U(zero  pressure)  . = A_  - Iu  - D/2  - B - e. . (22) 

mixture  P H 12  ' 

This  is  consistent  with  BM's  Figure  1 and  with  Stevenson's  Figure  1. 

(As  above,  the  B and  ei2  terms  may  be  neglected.)  Stevenson  obtains 
the  mixture  curve  from 

,P 

U(P)  = U(0)  - / P' (dV(P')/dP')dP'  . (23) 

o 

*D.  J.  Stevenson,  Nature  258,  222  (1975). 

17Ground- state  energy  (with  exchange)  plus  ep. 
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Figure  3.  Total  energy  U vs  volume  (as  smoothed  by  spline  routine) 


(24) 


V (P)  = V(NH3,P)  + JV(H2,P) 

where  V(NH3,P)  and  Vf^iP)  are  the  molecular  volume-pressure  relations 

of  ammonia  and  hydrogen  solids,  respectively.  Stevenson  obtained 

V (NH3,P)  by  interpolating  between  the  experimental  equation  of  state 

at  low  pressure  (P  < 20  kbar)  and  the  theoretical  equation  of  state 

at  high  pressure  (P  > 2 Mbar)  calculated1^  from  Thomas-Fermi-Dirac 

20 

theory.  V(H2,P)  is  taken  from  Ross's  analysis  of  the  Livermore 
shock  data. 

There  are  three  differences  between  our  metallic  NH4  U vs  V 
curve  and  BM's  curve: 

1.  Our  equilibrium  point  is  shifted  considerably  toward  Stevenson's 
value  of  - 95  x 10-2^  cm3/N  atom  (5.35  a.u.  radius).  We  obtain  - 76.7 
(4.99  a.u.)  whereas  BM  obtained21  47  (4.23  a.u.).  At  small  volume  (high 
pressure)  our  curve  tends  to  agree  with  BM's;  at  very  large  volume 

(low  pressure)  our  curve  approaches  Stevenson's. 

2.  At  the  respective  equilibrium  volumes  our  internal  energy  is 
about  0.5  ev.  higher  than  BM's.  BM  obtained  -5.84  ev.;  we  obtain 
-5.35  ev.  (Stevenson  has  -5.36  ev.).  The  0.5  ev  difference  between 
our  curve  and  BM's  exists  over  a sizable  part  of  the  volume  range. 

18J.  W.  Stewart,  J.  Chem.  Phys.  33,  128  (1956). 

*®E.  E.  Salpeter  and  H.  S.  Zapolsky,  Phys.  Rev.  158,  876  (1967). 

20M.  Ross,  J.  Chem.  Phys.  60,  3634  (1974). 

21Given  in  their  text;  their  Figure  1 indicates  approximately  52  (4.38  a.u 
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3.  Over  a considerable  part  of  the  range  of  interest  our  curve 
has  a smaller  slope  leading  to  smaller  pressure  at  a given  volume. 

We  agree  with  Stevenson  (Ref.  1)  that  BM's  transition  pressure 
(Pt)  estimate  is  invalid  due  to  their  not  obtaining  an  NH^  - j 
equation  of  state.  BM,  using  only  the  mixture  equilibrium  point , 

estimate  Pt  to  be  100  kbar  with  a range  from  60  to  140  kbar  based  on 

22 

an  0.5  ev.  uncertainty  in  A . 

In  computing  the  Gibbs  energy  (G  = U + PV  - TS;  T taken  as  zero), 
the  differences  (#2  and  #3  as  listed  above)  between  our  U vs  V and 
BM's  tend  to  cancel  and  we  obtain  essentially  the  same  G(P)  curve  as 
they  do.  At  very  low  pressure  our  curve  tends  to  approach  Stevenson's 
(see  Figure  4). 


VI.  SUMMARY  AND  CONCLUSIONS 

This  technical  report  is  primarily  concerned  with  the  BM 
calculation  for  metallic  NH^.  Our  calculation,  being  essentially  a 
repeat  of  BM's,  suffers  from  the  same  deficiencies,  namely: 

1.  The  calculation  is  essentially  atomic.  (Probably  minor.) 

2.  Sphericized  NH4+  ion.  (Probably  minor.) 

3.  Rigid  r0  (radius  of  hydrogen  "shell"). 

*D.  J.  Stevenson,  Nature  258,  222  (1975). 
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Our  total  energy  curve  (U  vs  V)  for  the  metal  would  lead  to  Pt  - 135 
kbar  using  the  BM  procedure;  we  emphasize  again  that  this  procedure 
is  invalid. 
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4.  Rigid  core  states  (Is,  2s,  2p) , merely  renormalized  for 
different  values  of  "a". 

5.  Approximation  made  in  using  BM's  iterative  procedure  for 
solving  their  eq.  (I)  with  exchange  (our  eq.  (12)). 

6.  Free  electron  calculation  for  the  Fermi  level.  (Probably 
significant,  particularly  at  the  higher  pressures.) 

The  major  deficiency  in  the  BM  paper  is  associated  with  their  attempt 
to  calculate  the  transition  pressure;  they  do  not  have  a curve  for  the 
NHj  - y H2  mixture! 

As  outlined  in  Section  IV,  our  attempts  to  duplicate  BM's  curve  I 
(ground-state  metallic  energy  without  exchange)  were  unsuccessful. 

Our  attempts  included:  Simply  solving  BM's  eq.  (I)  (our  eq.  (12))  with 
exchange  ignored  and  with  Xls  = X2s  = 0 (Sec.  IV.A.);  orthogonal! zing 
the  resultant  u(r)  to  u^s  and  u2s  and  then  finding  an  "expectation" 
energy  (Sec.  IV. B.);  including  X^s  and  X2s  using  a Rosenbrock  technique 
(Sec.  IV.C.).  Comparison  between  the  latter  calculation  and  BM's  is 
summarized  in  columns  three  and  four  of  Table  4 and  in  the  upper  part 
of  Figure  1. 

Attempts  to  match  BM's  calculation  including  exchange  were  more 
23 

successful  ; results  are  summarized  in  the  last  two  columns  of 
Table  4 and  in  the  lower  part  of  Figure  1.  We  emphasize  that  the 
validity  of  BM's  mathematical  approach  (solving  our  eq.  (12)  with 

2-*We  have,  however,  not  really  duplicated  BM's  results;  the  differences 
may  be  primarily  a reflection  of  the  approximately  two  decades  of 
development  in  electronic  computers  which  took  place  between  the 
two  calculations. 
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Ajj  = A2s  = 0 and  then  orthogonal izing  in  combination  with  an  iterative 
process)  has  not  really  been  established. 

Total  internal  energy  U is  given  by  the  sum  of  ground-state 
energy  (with  exchange)  and  mean  Fermi  energy.  Our  U(V)  curves  as 
well  as  BM's  and  Stevenson's  are  shown  in  Figure  3.  We  find  the 
internal  energy  at  the  equilibrium  volume  to  be  about  0.5  ev  higher 
than  that  obtained  by  BM,  i.e.,  ours  is  very  close  to  Stevenson's 
value  of  -5.36  ev.  We  note  that  Stevenson  uses  a proton  affinity 
value  some  0.46  ev  smaller  than  BM's';  this  would  tend  to  make 
Stevenson's  U(V)  curve  correspondingly  higher  than  BM's.  We  used  the 
BM  proton  affinity  value.  Stevenson's  equilibrium  radius  (for  the  metal) 
is  5.35  a.u. ; BM's  is  4.23  a.u.  We  obtain  4.99  a.u.,  i.e.,  shifted 
outward  considerably  toward  Stevenson's  value.  The  approximately 
0.5  ev  difference  between  our  metallic  U(V)  and  BM's  exists  over  a 
sizable  part  of  the  volume  range.  Our  U(V)  has  a smaller  slope  than 
BM's  over  a sizable  part  of  the  volume  range;  this  leads  to  smaller 
pressure  at  a given  volume. 

In  computing  the  Gibbs  energy  the  last  two  differences  tend  to 

cancel  and  we  obtain  essentially  the  same  G(P)  curve  as  BM.  At  very 

low  pressure  our  curve  tends  to  approach  Stevenson's  (see  Figure  4). 

1 4 

Stevenson  * seems  to  feel  fairly  strongly  that  the  desired 

4 

transition  will  not  take  place  and  has  suggested  that  BM  may  have 

*D.  J.  Stevenson,  Nature  258,  222  (1975). 

4 

D.  J.  Stevenson  (private  communication)  (Letter  to  Curtis  Selph, 

Oct.  1976). 
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made  a computational  error.  Our  repeat  of  the  BM  calculation  indicates 
one  cannot  explain  the  difference  between  Stevenson's  NH4  curve  and 
BM's  (see  Figure  4)  solely  by  computational  error  on  the  part  of  BM.^4 
As  an  argument  in  support  of  his  results  over  BM's,  Stevenson  cites4 
the  fact  that  his  equilibrium  volume  is  much  larger  than  BM's;  our 
calculation  shows  that,  within  the  BM  approach,  one  can  shift  the 
equilibrium  volume  considerably  without  making  an  appreciable  change 
in  G(P).  It  is  extremely  difficult  to  estimate  ''error  bars"  for  the 
curves  in  Figure  4 (or  Figure  3);  they  may  be  as  large  as  1-3  ev.  In 
our  opinion,  while  a transition  from  the  NHj-%H2  mixture  to  metallic 
NH4  in  the  funny  sodium  form  is  possible,  it  appears  unlikely,  par- 
ticularly below  about  one  Mbar.* 


4D.  J.  Stevenson  (private  communication)  (Letter  to  Curtis  Selph, 

Oct.  1976). 

^4This  assumes  that  the  BM  Is,  2s,  and  2p  wave  functions  are  essentially 
correct;  we  did  not  recalculate  these.  We  note  that  Stevenson  did 
his  calculation  twice;  once  with  the  BM  cores,  once  with  another  set  - 
the  difference  in  equilibrium  U was  only  0.04  ev.;  difference  in 
equilibrium  radius  only  0.17  a.u.  Thus,  recalculating  the  core  wave 
functions  is  unlikely  to  change  our  results. 

*Note  added  in  proof:  Discussion  with  Dr.  Stevenson  indicates  that 
funny  sodium  is  just  not  a good  candidate  for  the  metallic  form  of 
NH4  (assuming  his  mixture  curve  to  be  reasonably  good,  the  volume 
would  have  to  change  drastically  to  get  a common  U vs  V tangent).  It 
is  more  probable  to  go  to  an  arrangement  such  as  H's  forming  an  "fee" 
lattice  with  N sitting  in  the  "body-center"  position  or  possibly  to  a 
metallic  state  formed  by  the  overlapping  in  k- space  of  the  molecular 
energy  bands.  Stevenson  estimates  the  required  pressures  to  be  in  the 
Mbar  (and  above)  range.  After  discussion  of  our  results  (see,  in  par- 
ticular, the  last  paragraph  of  this  report)  Stevenson  indicated  that 
BM  may  have  an  error  in  their  physics  (e.g.,  use  of  the  free-electron 
Ep  can  lead  to  significant  error) . 
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APPENDIX  A 


DERIVATION  OF  THE  ONE-VARIABLE-OPERATOR 
INTEGRALS  OF  THE  HARTREE-FOCK  EQUATION 


Explicitly,  the  N-electron  wave  function  ¥ of  the  main  text  is 


^(x^)  • * • ^j(xN) 

*2c*i)  ' ’ ‘ ^2^ 


(A-l) 


<t>N  (xi) 


with  ^(Xj)  = 
si(?j)  can  be 


. <j>^(r\)  is  a purely  spatial  function, 

either  a(6j)  or  B(5j)  where 


a(Oj)  = 1,  if  the  spin- coordinate  a-  is  up. 
=0,  if  Oj  is  down. 

B(Oj)  = 0,  if  is  up. 

= 1,  if  Oj  is  down. 


The  orthogonality  condition  (eq.  (5)  of  the  main  text)  written 
explicitly  for  space  and  spin  is 


C1!.’ dT  hSjji  s k CoSl  > sm  tCF*,’  ■ *m 


(A-2) 


6^  is  a condition  to  be  imposed  on  the  spatial  functions;  6S , sm  is 
insured  by  the  a, 6 definitions  given  above. 

Consider  Hjc^net^c  of  eq.  (4)  of  the  main  text;  this  one-variable 
operator  acts  on  det  [^(r^)].  (Since  we  deal  here  with  a one-variable 
operator,  and  since  the  operator  does  not  affect  the  spin,  we  can 
ignore  the  spin  coordinate;  in  Appendix  B,  where  we  deal  with  a two- 
variable  operator,  we  explicitly  include  the  spin  functions.)  Now 
consider  just  one  term  of  Hjc^net^c,  namely  V*:  ignoring  the  purely 
multiplicative  factor  this  operator,  acting  on  det  Crj ) ] gives 
N ^ 

l V*<J>  (rj  [permutations  of  <{>  (£)]  . (A-3) 

j=l  1 3 1 a b J 


For  each  value  of  j the  square  bracket  contains  (N-l)f  terms,  each  with 
(N-l)  factors:  a = 1,2, ’“.N  (a  ^ j);  b = 2,3,***,N.  Combining  this 
with  det*[4>^(rjc)]  and  using  the  orthogonality  conditions  we  see  that 
only  one  of  the  N!  terms  of  det*  "matches"  (gives  a non-zero  integral) 
any  particular  term  of  (A-3),  i.e., 

4>j  (r^)  [correct  permutation  of  4>*(^b)]j 


matches 


VJ* . (?i)  [particular  permutation  of  4»a ] j 


This  leads  (for  V*  alone)  to 


(N 


-1)1  I /♦.*(r1)VV.(?,)dT1  • 

j=l  1 i l J i i 


(A-4) 
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We  obtain  similar  expressions  for  V*  V2,,,,t  V2.  Since  the  variable 

2 3 N 

acted  on  in  (A-4)  is  the  variable  of  integration,  each  of  V2,  V2,  etc., 
will  give  the  same  result.  Thus  (dividing  by  the  N!  of  the  denominator): 

^kineticl^/^l^  ■ -C«2/2m)Ji  J<j>  * (?)Vj<J> . (A-5) 

This  is  (A)  of  the  main  text. 


A similar  treatment  of  H 

nuclear 


gives  (B)  of  the  main  text. 


APPENDIX  B 


DERIVATION  OF  THE  TWO-VARIABLE-OPERATOR 
INTEGRALS  OF  THE  HARTREE-FOCK  EQUATION 


The  third  term  of  eq.  (4)  of  the  main  text,  the  electron-electron 
interaction  He£_e£>  is  a two-variable  operator.  Consider  just  one  term 
of  this  operator,  namely  l/lr^-r^l  , operating  on  det  [^(x^)].  This 
gives 


N 

l l 

k=l  £=i 


I ^2  I 


sk(^i)s£(02)  [Permutations  of  <t'a(r|,)sa(aj))  ] 


k£ 

(B—  1 ) 


For  each  choice  of  k,£  in  the  double  sum,  the  square  bracket  contains 
(N-2)I  terms,  each  with  (N-2)  factors:  a = 1,2,3,***,N  (a  j*  k,£); 
b = 3,4,*-,,N.  when  we  multiply  by  det*[<j>^(xj)] , two  of  the  N!  terms 
in  det*  will  "match"  any  particular  term  of  (B-l),  i.e., 

<l>^(r1)^(r2)sk(^1)sJl(?f2)[correct  permutation  of  <J>* (rb)sa(CT^) ]kj^ 
and 


M*l^k^H^l)sk(^tcorrect  Permutation  of  4>a(^)sa(ab) ] 
match 

— |'J.  ITj sJc(a1)s£(o2)  [particular  permutation  of  4>a Cr*bD sg  (CTb)  ] 

as  far  as  spatial  functions  are  concerned.  For  the  spin  functions: 
s^,  s£  can  be  any  a, 3 combination  for  the  first  of  the  two  matching 
terms;  in  the  second  matching  term  s£  must  equal  s^  in  order  to  have  a 
non-zero  value;  this  is  readily  seen  to  imply  that,  for  this  term, 
electrons  1 and  2 must  have  parallel  spins.  Thus,  l/|i^-^|  gives 


1 

i 

« 

; 

i 
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w-wlUl  IV&TW 

« l?rf2l  1 ‘ 


- 6 (spins)  / 


I rl_r2 1 


dTldT2^‘ 


So  far,  we've  considered  only  l/|i^-r^|;  there  will  be  a total  of 
N(N-l)  operators: 

i/ifj-'2i.  i/ifj-fji,---,  i/ifrjji,  i/ir2-^i,— , i/i?N-fN.1i. 

r^,  rj  in  each  l/|r^-x^|  are  also  the  variables  of  integration  in  that 
particular  case;  thus,  each  operator  gives  the  same  result.  Then 
(dividing  by  the  N!  of  the  denominator  and  replacing  the  ez/2 
multiplicative  factor)  we  obtain  (C)  of  the  main  text.  (We've 
interchanged  r^,  r^  in  the  second  integral  of  (B-2);  this  is  clearly 
allowable  due  to  the  form  of  the  operator.) 
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APPENDIX  C 


CAN  THE  A MATRIX  APPEARING  IN  THE 
FOCK  EQUATIONS  ALWAYS  BE  DIAGONALIZED? 

There  is  considerable  confusion  on  this  point  in  the  literature. 

25 

Reitz  states  that  one  can  always  choose  solutions  such  that  A is  a 

diagonal  matrix  (emphasis  added).  Reitz  is  dealing  explicitly  with  the 

crystal  case.  Slater,  in  his  derivation  of  the  Fock  equations  for  the 
26 

atomic  case  , appears  to  show  that  a unitary  transformation  of  the 

one-electron  functions  <J>j  (q^)  can  be  made  which  will  diagonalize  A 

(Section  17-1);  however,  in  Section  17-5  (for  non-closed  shells)  of 

ref.  26  Slater  states  that  one  cannot  diagonalize  A.  In  the  crystal 
27 

case  Slater  seems  to  indicate,  in  Section  1-2,  that  one  can  diagonalize 
A. 

Without  pursuing  the  matter  further  it  would  appear  that  the 
crucial  factor  is  whether  one  uses  a single  determinantal  function 
(A  diagonalizable)  or  whether  one  needs  more  than  one  determinantal 
function  (A  not,  in  general,  diagonalizable). 


R.  Reitz  in  "Solid  State  Physics,"  ed.  Seitz  and  Turnbull,  Vol.  1, 
p.  1,  Academic  Press  (1955). 

2<\J.  C.  Slater,  "Quantum  Theory  of  Atomic  Structure,"  Vol.  2,  McGraw- 
Hill  (1960). 

77J.  C.  Slater,  "Quantum  Theory  of  Molecules  and  Solids,"  Vol.  4, 
McGraw-Hill  (1974). 
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APPENDIX  D 


DERIVATION  OF  THE  BERNAL-MASSEY  EQUATION: 
FIRST  STEPS  IN  OBTAINING  THE  RADIAL  EQUATION 


In  spherical  coordinates 

V2  = (l/r2)3/3r(r23/3r)  + (1/r2) [ (l/sin0) 
r 


3/3e(sin63/36)  + (l/sin26)32/3<t>2] 


We  knovr0  that  the  square  bracket  of  (D-l)  operating  on  Y£m  gives 
-£.(£+l)Y_  . Since  H = 0 here,  we  may  drop  this  term. 


Eq.  (8)  of  the  main  text  then  becomes 

. *2  Y0  d2u(ri)  Y0  ^ #Wrl>  ~ 

2m  "dr/'  V°(rl)  ^ u(rl}  + 6 u 


Cri)Y0 


ls,2s,2p  (r2) u(r2)Y0(92,<|)2) 

-e2  I 6 (spins) [/  J A 2A  ° dTgfy  (r.)  CD-2) 

j |rrr2|r2 


u(ri) 


Y0  " Xls 


Uls(rl) 


Y0  " X2s 


U2s(rl) 


Since  the  YQ  terms  have  no  angular  dependence  and  as  all  are  normalized 
with  the  same  factor  l//frr,  we  may  simply  factor  out  the  YQ;  factoring 
out  1/r^  as  well  and  substituting  for  <j>j  as  per  eq.  (9)  of  the  main 
text,  we  obtain  eq.  (10)  of  the  main  text.  The  direct  coulomb  term 
has  been  replaced  by  an  averaged  term  Ze£f(rj)/rj. 


Schiff,  "Quantum  Mechanics,"  McGraw-Hill  (1949).  See  eqs.  (14.21) 
and  (14.22),  Chap.  IV. 
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APPENDIX  E 


DERIVATION  OF  THE  BERNAL-MASSEY  EQUATION: 
TRANSFORMING  TO  ATOMIC  UNITS 


Let 

r = a0p 

with 

a0  = “fi/me2  = 1 Bohr  radius 
,fi2/2ma02  = 1 Rydberg 

Also  let 

VoCPj)  = VoOj) 
uCP)  = u(r) 
zeff(Pl)  = Zeff^ri) 

After  all  these  substitutions  have  been  made  and  the  entire  equation 
multiplied  by  2ma02/-ft2,  replace  pj  and  P2  with  r and  rj  respectively. 
This  gives  eq.  (11)  of  the  main  text. 
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; 

; 

li! 

APPENDIX  F 

DERIVATION  OF  THE  BERNAL-MASSEY  EQUATION: 

i 

MANIPULATION  OF  THE  EXCHANGE  TERM  FOR  CLOSED  £-SUBSHELLS 

Due  to  the  rotational  invariance  in  eq.  (11)  of  the  main  text  we 
can  always  (in  the  integration  over  r^)  choose  the  polar  axis  to  be  in 
the  r direction. 

t'  = (r2  + rj2  - 2rrj  cosy)*s.  (F-l) 
For  rj  < r, 

r'  = r[l  + (rj/r)2  - 2(r1/r)cosy]Js 
= rX, 

OO 

with  X"1  = l PjT(cosy)  (r./r)£  . 

£'=  0 

(ref.  28,  eq.  (14.10))  (F-2) 

Similarly,  for  r^  > r, 

00  > 

I P.-fcosyJfr/rj)1  . (F-3) 

r rir=o  1 1 

For  j standing  for  Is,  the  r^  < r part  of  the  exchange  integral  in  eq. 
(11)  of  the  main  text  is  then 

1 r 00  o - / 2ir 

[7  / uis(ri)u(ri)  I (*i/r)  drx  / P ^ (cosy) d (cosy)  — / d^,]  (F-4) 

r r^O  15  l'=0  y=0  Jf*  <p1=0  1 

using  Yls  = 1//Itt  . Inserting  Pq(cosy)  = 1 into  (F-4)  and  using 

(ref.  28,  eq.  (14.15)), 

^®L.  Schiff,  "Quantum  Mechanics,"  McGraw-Hill  (1949). 
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CF-S) 


ir 

/ P.  (cosy) P„' (cosy) d (cosy)  = (2/[2£+l])6t  , 

Y=0  x 


the  Is  exchange  term  (for  r^  < r)  in  eq.  (11)  becomes 


■ (2/r) [ /r^uls(r1)u(r1)dr1]uls(r)  . (F-6) 

For  rj  > r we  have 

-2[  / uls(r1)u(r1)(l/r1)dr1]uls(r)  . (F-7) 


There  will  be  similar  expressions  for  the  2s  term. 
For  j standing  for  2p 


£ 

l Ylm(6l^l)Ylm(6*4,)  = (3/47T)  cos9l  cos0 
m=-£ 

+ (3/8ir)sin0^  sin0[e2  ^”^1^  + e"^^-(^l^] 


= (3/4it)[cos0^  cos0  ♦ sinOj  sin0  cos (<j>-<J>^) ] 

* (3/4tt)cosy  = (3/4tt)P1  (cosy)  . (F-8) 

The  2p  integral  in  eq.  (11)  is  then  (for  r^  < r)  , 

r “ tt 

[l/r  / u2n(Ti)u(Tl)  I (Ti/T)  dri  / (3/4n)P. (COSY) 
rj-0  P £'-0  Y=0 

2ir 

X P0>  (cosy)  d (cosy)  / d<J>  ] . 

* *r=° 


Using  (F-5) , the  2p  exchange  term  (for  r < r)  becomes 


C-2/r2)[  / u2p(r1)u(r1)r1dr1]u2p(r)  , (F-9) 

r,  =0 
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Collecting  eqs.  (F-6)  through  (F-10),  the  entire  exchange  term  can  be 


written 


-2  l {r"*"1/  u^O^Mr  Jr^drj 
ls,2s,2p  o 

o 3.  -Z- 1 

+ r / V(rl)u(rl)rl"  " drl}un*(r)  ' 

r 


CF— 11) 


Inserting  this  into  eq.  (11)  of  the  main  text  and  rearranging  somewhat 
we  obtain  eq.  (12)  of  the  main  text. 


APPENDIX  G 


THE  NUMEROV  TECHNIQUE  AND  COMPUTER  PROGRAM  "WAVEB" 

I.  Basis  of  the  Numerov  Process  and  a Narrative  Description  of  "WAVEB" 
We  note  that  "MAIN"  of  the  computer  program  "COMBINE"  is  very 
similar  to  "WAVEB";  the  differences  are  discussed  in  Part  IV  of  this 
Appendix. 
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The  treatment  here  is  based  on  that  of  Hartree  . The  Numerov 
process  is  applicable  to  linear  second-order  differential  equations. 
Consider  the  equation 

u"  = f (r)u  * g(r)  (G-l) 

We  treat  three  cases: 

(1)  No  exchange;  Xjs  = X£s  = 0.  Eq.  (13)  of  main  text. 

(2)  No  exchange;  X's  included.  Eq.  (16)  of  main  text. 

(3)  With  exchange;  X^s  = X2S  = Eq.  (18)  of  main  text. 


In  all  three 

cases: 

f(r)  = V(r)  - e . 

For  case 

(1), 

g(r)  = 0 

(G-la) 

For  case 

(2), 

g(r)  = XlsUls(r)  + X2su2s(r) 

(G-lb) 

For  case 

(3), 

g(r)  = l (-2){r-X'_1/ru  „(x)u  (x)x£dx 

ls,2s,2p  o 

♦ rlf  unJl(x)uw(x)x'*'‘1dx}unJl(r) 

(G-lc) 

^®D.  R.  Hartree,  "The  Calculation  of  Atomic  Structures,"  Wiley  (1957), 
Chapter  4. 
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where  u^(x)  is  the  u of  the  previous  go-round  when  we  are  solving 
eq.  (G-l)  for  U(k+1)(r). 

Define  the  "first  difference"  as 


6u-  i = u . , - u. 

J+5s  j+l  j 


(G-2) 


where  Uj  is  a shorthand  notation  for  u(r^)  where  rj  here  represents 
the  value  of  r at  the  jth  point  of  the  r mesh.  Define  the  "second 
difference"  as 


= 5 Vis  " 6uj-*s  = “j+1 

A Taylor  series  expansion  of  u leads  to 


- 2u.  + u. 


j-l 


(G-3) 


62Uj  = (Ar) 2[Uj"  ♦ ± 52Uj"  - ^ 64Uj"]  + O(Ar)8  . (G-4) 

We  will  ignore  the  A^term  and  the  O(Ar)8  term  in  eq.  (G-4).  Combining 
eqs.  (G-l),  (G-3),  and  (G-4)  gives 


Vi  = rTf1'  ~T  {uj[2S+10fj3  - uj-ils-€j-0  + 8j+i  + 10«j  + gj-i} 

1 j+r 

(G-5) 

with  S = 12/ (Ar) 2 

An  inspection  of  eq.  (G-5)  shows  that  u^+1  depends  on  the  u values 
of  the  two  "previous"  (further  left)  points  j and  j-l.  Thus,  some  other 
procedure  must  be  used  to  determine  Uj  near  r = 0.  This  "starting" 
procedure  is  based  on  the  expansion 


un^(r)  = Arx'+1(l  + ar  + gr2  + yr3  + *••) 


(G-6) 


near  r = 0.  Substitution  of  (G-6)  into  (G-l)  and  solution  in  series 
allows  evaluation  of  a,  g,  and  y in  terms  of  f and  g.  This  is  used  to 
generate  u(2)  and  u(3).  (From  (G-6),  u(0)  = 0.0.)  The  explicit 
evaluation  of  a,  g,  y is  given  in  Part  II  of  this  appendix. 
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WAVEB  may  be  thought  of  as  a two-part  process:  the  first  involves 
getting  the  correct  number  of  nodes;  the  second  involves  getting  zero 
slope  at  the  r = a end  of  the  mesh.  Both  parts  involve  a number  of 
iterations  (typically  on  the  order  of  30  altogether). 

Part  one  of  WAVEB  proceeds  by  an  iterative  process  as  follows: 

An  initial  guess  is  made  for  e (read  in  as  input  data).  u(2)  and  u(3) 
are  then  found  using  eq.  (G-6).  The  remaining  u(r)  are  then  generated 
(for  j = 4 to  j = N)  using  eq.  (G-5) . An  atomic-like  function  labeled 
by  quantum  numbers  n and  i must  have  n-5,-1  nodes.  Once  u(r)  for  j = 

1,N  is  generated,  the  number  of  nodes  are  counted  and  e is  changed  to 
e+Ae  or  to  e-Ae  to  start  the  next  iteration.  From  the  form  of  eq. 

(G-l)  and  either  eq.  (13)  or  (16)  of  the  main  text  one  sees  that  e 
should  be  made  more  positive  if  there  are  too  few  nodes  (need  more 
curvature)  etc.  The  magnitude  of  the  initial  Ae  step  is  read  in  as 
data;  as  subsequent  iterations  "zero-in"  to  just  produce  the  right 
number  of  nodes  within  the  prescribed  mesh,  the  size  of  Ae  is  reduced. 

• When  the  correct  number  of  nodes  are  found  (with  the  last  node 
sufficiently  close  to  the  r = a point)  the  program  switches  to  part 
two  - the  slope  test.  The  form  of  eqs.  (G-l)  and  (13)  or  (16)  indicate 
that  when 

(slope) (sign  of  u(N)) 

is  negative  (indicates  slope  is  toward  the  horizontal  axis;  less 
curvature  is  required)  one  should  decrease  e;  when  positive,  one  should 
increase  e.  This  process  is  continued  until  | Ac | becomes  smaller  than 
some  preset  value. 
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11 • Evaluation  of  a,  g,  and  y for  the  "Starting"  Procedure 
For  u3s,  eq.  (G-6)  is  (for  r near  r = 0) , 

u(k+l)(r)  = A(k+l)(r  + ar2  + 6r3  + Yr^  ' (G-7) 

From  eq.  (12)  of  the  main  text, 

V(r)  = V0(r)  + 2Zeff(r)/r 

= -14/r  - 8/r0  + 2Ze££(r)/r,  (r  near  r = 0) . 

For  the  three  cases  listed  in  connection  with  eq.  (G-l) 

f(r)  = -14/r  - 8/r0  + 2Zef£(r)/r  - e . (G-8) 

Expanding  Zeff(r)  in  a taylor  series  (for  r near  r = 0) , 

f (r)  = -14/r  + 2Zef£(0)  - 8/r0  - e + rZgf£(0)  (G-9) 

f H 

where  we  have  used  Z(0)  = 0.0.  Ze££(0)  and  Zg££(0),  as  well  as 

and  Ansans  introduced  later,  may  be  evaluated  by  expanding  the 
appropriate  function  (Ze££(r)  or  u^OO)  in  a Taylor  series  for  the 
first  three  non-zero  values  of  the  function  and  solving  the  resultant 
equations  simultaneously.  We  use  the  Zg££(r)  of  BM's  Table  I.  We  g 

* tf 

obtain  Zg££ (0)  = 20.1340  and  Zg££(0)  = -17.4959;  rQ  = 1.84  as  per  BM. 

Substituting  (G-7)  and  (G-9)  into  (G-l)  we  obtain 
A(k+i)C2a  + 6&r  + 12yr2]  - g(r) 

+ [14/r  + R£  - rZeff(0)]A(k+1)[r  + ar2  + gr3  + yr1*]  = 0 (G-10) 

with 

R£  = e + 8/r0  - 2Zgff(0)  (G-ll) 
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Case  (1):  g(r)  = 0 


Equating  coefficients  of  the  various  powers  of  r to  zero  separately 
(keeping  terms  to  order  r2  only) , we  obtain 


a = -7 


6 = -[14a  + Re]/6 
y = [Zeff(0)  - 14g  - <xRe]/12 


(G- 12) 


Cases  (2)  and  (3) 
From  (G-6), 


uis(r)  - Als(r  + alsr2  ♦ blsr3  + c^r") 

u2s(r)  = A2s(r  + a2sr2  + b2sr3  + c2sr^ 

u2pW  = A2p(r2  + a2pr3  + b2pr4  + c2pr^ 


To  order  r , g(r)  can  be  written 


(G-13) 


g(r)  = AlsMls(r  + alsr2)  + A2sM2s(r  + a2sr2)  . (G-14) 
For  case  (2),  is  clearly  just  Ans.  In  Part  III  of  this  appendix 
we  show  that  g(r)  for  case  (3)  can  be  written  as  in  eq.  (G-14)  with 


Mns  = -2/  unsWu(k)Wx'ldx  * 
o v 


(G-15) 


for  the  u^+i)  solution.  for  case  (3)  is  evaluated  in  Subroutine 

"AUX2"  of  "COMBINE".  (We  note  that  the  U2p  function  does  not  appear 
in  f(r)  for  r near  r = 0.)  When  (G-14)  is  substituted  into  (G-10)  and 
the  resultant  equation  is  reordered  we  obtain 


\k+i)E2a  ♦ 14J  + {Vi)Cl4a  + 6B  + Re3  - l C'WW** 


(G-16) 


♦ {A(k+l)taRe  + 14B  + 12^  - zeff(°)]  ‘ l CAnsMnsans)>r2  = 0 
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Equating  coefficients  of  the  various  powers  of  r to  zero  separately 
(keeping  terms  to  order  r2  only) , we  obtain 


a = -7 

2 

B = -[14c  * Re  - A'*  (G-17) 

v n=l 

1 ' - aRe  - 146  * l (\s<\s^s)Vl2  • 

v n=l 

III.  Evaluation  of  for  Case  (3). 

Using  eq.  (G-7)  for  u^  (x)  and  eq.  (G- 13)  for  u^,  the  first 
integral  in  eq.  (G-lc)  becomes 

^£A(k)[l  + an£r  + Vr'  + 

T 

X / xJ,+1[l  + a^x  ♦ bn£x2  + ***]x[l  + ax  + 8x2  + "’Ix^dx  . 
o 

(The  r"A-l  and  un^(r)  factors  of  (G-lc)  have  been  included  in  the 
expression  above.)  The  lowest-power  term  is 

**„/>**  - . tG-18) 

We  are  interested  only  in  terms  of  0(r2)  or  lower  (see  eq.  (G-16)). 
Thus,  even  for  l = 0,  this  first  integral  term  may  be  ignored. 

We  rewrite  the  second  integral  term  in  eq.  (G-lc)  as 

An£(r)[/a  - n (G-19) 

o o 

The  first  integral  in  (G-19)  must  be  a constant  (independent  of  r) ; 
writing  for  this  constant,  this  first  term  of  (G-19)  becomes 

T\zT^l  + an£r  + "‘Kl 
using  eq.  (G- 13) ; we  rearrange  this  as 
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(G-20) 


An£Cn£,r 


2SL+1 


[1 


+ an£,r  + 


For  U2p(Jl=l),  this  will  be  0(r3)  and  higher:  thus  U2p  can  be  ignored 
here.  For  i = 0 we  have 


AC  r + ACar2  + 
ns  ns  ns  ns  ns 


(G-21) 


which  must  be  considered.  Using  eq.  (G-7)  for  and  eq.  (G- 13) 

for  u^,  the  second  integral  in  (G- 19)  becomes 

+ an4r  + *'•] 


X / x£+1[l  + a^x  + •••]x[l  + ax  + •••]x  * *dx  . 


-1-1 


o 

0 

(The  r u^fr)  of  eq.  (G-lc)  or  (G- 19)  is  included  in  this  expression.) 
The  lowest-power  term  is 


-AZzA(k)r2Z+1/*xdx  = -An£A(k)r2,l+3/2  (G-22^ 

which  is  of  0 (r3)  even  for  X.  = 0;  thus,  this  term  may  be  ignored. 

Thus,  contributions  to  the  "starting"  equation  from  g(r)  come 
only  from  u^s  and  U2S  as  indicated  by  (G-21);  this  gives,  for  r near 
r = 0, 

2 a i 

g(r)  = -2  l A ns[/  unS(x)urk) (x)x_1dx](r  + ansr2)  (G-23) 

n=l  o 

verifying  eqs.  (G-14)  and  (G-15)  for  case  (3),  i.e.,  for  the  exchange 
case  with  X^s  = X2S  = 0. 
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IV.  Difference  Between  MAIN  of  "COMBINE"  and  "WAVEB" 

The  overall  routing  of  "COMBINE"  is  shown  in  Figure  5.  The 
discussion  in  Parts  I,  II  and  III  of  this  appendix  applies  to  MAIN  of 
COMBINE  as  well  as  to  WAVEB. 

Our  older  system  for  generating  an  atomic-like  wave  function 
worked  reasonably  well  on  a single-shot  basis  but  isn't  good  for  the  BM 
problem  where  we  "recycle",  i.e.,  where  the  whole  iteration  is  repeated 
again  with  a new  g(r).  (In  this  discussion  the  word  "iteration" 
refers  to  a new  guess  at  the  energy  e for  a given  g(r);  the  word 
"cycle"  refers  to  the  set  of  e iterations  for  a fixed  g(r).) 

The  specific  problem  with  the  older  system  is  the  result  of  several 
factors.  By  the  third  or  fourth  cycle  "convergence"  is  fairly  well 
obtained  (input  u (r)  z output  u (r)).  Under  the  old  system,  u(r) 
is  forced  to  have  an  extra  node,  then  energy  is  decreased  somewhat  and 
e,u(r)  are  found  which  give  zero  slope.  After,  say,  the  third  cycle, 
the  input  energy  for  the  first  iteration  of  the  fourth  cycle  is  just 
about  right  to  give  zero  slope  with  the  correct  number  of  nodes; 
however,  the  old  system  forces  an  increase  in  energy  to  create  an 
extra  node,  i.e.,  forces  one  away  from  the  proper  u(r). 

In  the  new  system  the  program:  (1)  Gets  the  correct  number  of 
nodes.  (2)  Increases  the  energy,  pushing  the  last  (greatest  r value) 
node  inward  until  this  last  node  is  some  preset  number  of  spaces  in 
from  the  end  of  the  mesh.  (3)  Juggles  the  energy  to  give  zero  slope 
at  r = a.  Since  the  new  system  needs  only  25  - 35  iterations/cycle, 
we  start  each  cycle  with  the  original  | Ae | We  have  also  put  the 
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DE-DEMIN  test  (the  test  for  proper  exit  from  the  e-iteration  loop)  on 
both  the  positive  Ae  and  negative  Ae  sides. 


Results  of  these  changes:  the  first  real  run  of  COMBINE  using  the 
old  system  involved  five  cycles;  the  fifth  cycle  required  292  iterations! 
With  the  new  system  this  cycle  obtained  virtually  the  same  final  e 
(differed  by  1 in  the  7th  significant  digit)  using  only  29  iterations. 
Secondly,  using  a starting  Ae  for  this  fifth  cycle  equal  to  (original 
Ae)/10  only  shortened  the  number  of  iterations  by  6;  the  original  Ae 
(equal  to  0.1  Ry  in  this  case)  gave  virtually  the  same  final  e and 
final  Ae  and  produced  a small  enough  final  slope;  thus  the  new  system 
just  uses  the  same  starting  Ae  (of  order  of  0.1  Ry)  for  every  cycle. 
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THE  "TWAVE"  PROGRAM 

I.  A Narrative  Description  of  "TWAVE" 

The  heart  of  the  3-dimensional  search  is  in  Subroutine  DRSNIN,  a 
Ray  Scanlon  routine  based  on  the  Rosenbrock  technique.  Initial  guesses 
for  Xj,  X2,  and  £ and  the  initial  step  sizes  for  each  are  read  in  as 
data.  The  unnormalized  BM  uls(r)  and  u2s(r)  functions  are  read  in  and 
normalized. 

DRSNIN  contains  a sophisticated  search  procedure.  For  each 
"dimension"  the  routine  first  steps  in  the  positive  direction,  then  in 
the  negative  direction  if  necessary;  a success  (smaller  "object" 
function)  followed  by  a failure  (larger  object  function)  terminates 
the  search  in  that  direction.  Step  sizes  are  automatically  adjusted. 
Thus,  DRSNIN  attempts  to  minimize  the  object  function  which  in  our 
case  is 


(weight  factor) [du(r) /dr | a]2  + <U|UH>  “ * /U|U2S>  2 


(H—  1 ) 


This  reflects  the  boundary  conditions:  zero  slope  at  r = a;  u orthogonal 
to  Ujs  and  u2g.  The  object  function  is  evaluated  in  Function  (Subroutine) 
H(X) . H(X)  calls  Subroutine  LWAVE  which  uses  the  Numerov  technique 
(see  App.  G)  to  find  u(r).  When  u(r)  contains  the  wrong  number  of 
nodes  a large  number  (10000.)  is  substituted  for  du(r)/dr|a  in  the 
object  function.  The  program  can  be  run  in  either  of  two  modes: 

(1)  Terminates  after  a preset  number  of  cycles  or  when  the  object 
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function  becomes  smaller  than  some  preset  value,  which  ever  occurs 
first.  (2)  Terminates  after  a preset  number  of  cycles.  If  mode  (2)  is 
used,  one  gets  an  indication  of  the  convergence  from  a simple  inspection 
of  the  print-out.  It  should  be  noted  that,  typically,  the  object 
function  is  evaluated  many  times  within  one  "cycle". 


II.  Normalization  of  u(r) 

Because  of  the  presence  of  the  X u^ (r)  terms,  eq.  (16)  of  the 


main  text  is  not  a standard  eigenvalue  equation.  Defining 

6 = -d2/dr2  + V(r)  , 
g(r)  = X°u,  (r)  + X°u-  (r)  , 


we  may  write  eq.  (16)  as 


6u(r)  + g(r)  = eu(r) 


(H-2) 

(H-3) 


(H-4) 


Suppose  eq.  (H-4)  has  been  solved  to  obtain  e,  u(r).  Now  define 

U(r)  = u(r)/N  (H-5) 

with  N=  ^u^^  so  that  ^j|u^  =1.  Substituting  (H-5)  into  (H-4)  we 


obtain 


5u(r)  + g(r)/N  = eU(r)  . 


If  we  now  define  a new  G(r): 

G(r)  = X1uls(r)  + X2u2s(r)  (H-7) 

with  A-=A°/N,  so  that  G(r)=g(r)/N, 

flu(r)  + G(r)  = eU  (r)  (H-8) 

showing  that  U(r)  satisfies  essentially  the  same  equation  as  u(r), 
with  the  same  e,  but  with  new  values  of  A^  and  A2. 


t 
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COMMENTS  ON  THE  VALIDITY  OF  THE  BERNAL-MASSEY 
PROCEDURE  FOR  SOLUTION  OF  THEIR  EQUATION  (I) 

Within  the  approximations  of  the  BM  model  (i.e.,  solving  the  ion 
problem  first,  then  adding  4»-js C r)  = u(r)  while  keeping  the  core 
functions  fixed)  the  U2S(r)  equation  is 

-u2s(r)  + V2s(r)u2s(r)  ' e2su2s(r) 

+ Xls,2suls^  = -82s  W r-V 

where 

V2s(r)  " VoW  + 2Zeff(r)/r  ’ 

Zeff^  = 2ZlsM  + Z2s^  + 6Z2p^  » 
g2sCr)  = g(u2s'uls^  + 8(u2s,u2p^  ' 

with 

E -2^'<l'1/runS(x)un,(>(x)xtdx 
o 

+ un£Wun'o(x)x"*'"ldx}un£Cr)  ' 
r 

and  VQ(r)  is  the  potential  due  to  the  nucleus  plus  the  four  hydrogen 
protons.  Eq.  (I)  of  BM  (in  our  atomic  units)  is 

-u"(r)  + V3s(r)u(r)  - e3su(r)  ♦ *3s>lsuls(r) 

+ X3s,2su2s^  = -g3s(r)  ' (I)  of  BM 
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with 

V3sm  = Vo(r)  . 2z5|f(r)/r  , 
zef£(r)  = 2Zjg(r)  . 2Z2s(t)  * 6Z2p(T)  , 
g3s(r)  = g(u,uls)  + g(u,u2s)  + g(u,u2p)  . 

We  note  that 

V3s(r)  " V2s(r)  = 2Z2s(r)/r  * U'2) 

In  the  BM  iterative  procedure  for  determining  u(r)  and  e3s,  the 

"last  go-round"  equation  looks  like 

[-d2/dr2  + V3s(r)  - e3s]uU(r)  = -g3^(r)  (1-3) 

where  g (r)  uses  u -*•  (r)  instead  of  uu(r).  The  superscript  u means 
"unnormal ized";  the  superscript  X means  "orthogonal ized,  normalized". 
Explicitly: 

ux  (r)  = [uu (r)  + au2s(r)]/N  , (1-4) 

N = <^uu  + au2s|uu  + » 

“ = - (?U\u2s}  * 

The  significant  differences  to  note  between  BM's  eq.  (I)  and  our  eq. 
(1-3)  are  that  eq.  (1-3)  has  the  X's  set  equal  to  zero  and  has  an 
unorthogonalized  and  unnormalized  "solution  u"  (uu(r))  different  from 
the  u-4.  (r)  used  in  getting  the  exchange  term  g3x(r).  (We've 
assumed  "convergence",  i.e.,  u-^  (r)  of  the  jth  cycle  z u (r)  of 
the  (j-l)th  cycle.) 

We  adopt  the  point  of  view  that:  if  ux  (r)  is  a "proper  solution" 
to  the  problem  it  should  satisfy  eq.  (1)  of  BM;  we  first  rewrite  (I): 
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[-d2/dr2  + V3s(r)]u(r)  = e3j.u(r) 

-X3s,2sU2s(r)  ‘ g3sW  • V-V 

(We've  set  \3s  ^ = 0 as  did  BM;  our  numerical  work  also  indicates 
^lsluU^  *s  sma11*)  We  now  Put  u x (r)  into  the  left-hand-side  (LHS) 
of  eq.  (1-5)  and  see  if  we  can  obtain  the  RHS: 

LHS  of  (1-5)  « [-d2/dr2  + V3s(r)]ux  (r) 

= (1/N) [-d2uu/dr2  + V3suu] 

+ (a/N)[-d2u2s/dr2  + V3su2s]  , 

= (1/N)[e3suu  - g3f  ] 

+ (a/N)tCe2s+V3s-V2s)u2s-g2s-Xls,2suls^  » 
using  eqs.  (1-1)  and  (1-3).  With  some  manipulation  we  obtain 

LHS  - e3ju±  - (o./N)[e3s  - £2s]u2s  - g3f  (1-6) 


* C°/n)[(v3,  - v2s)„2s  - g2s  - qs,2suls]-[(l-N)/N]gA  . 
This  equals  the  RHS  of  (1-5)  if: 

1‘  X3s,2s  (a/N) [£3s  ‘ G2s^ 


2.  The  second  line  of  (1-6)  0. 


TABLE  1-1.  NUMERICAL  EVALUATION  OF  "EXTRA"  TERMS  IN  EQ.  (1-6) 
This  evaluation  is  for  a = 3.68. 
e3s  = "°-507»  a = 0.302,  N = 0.629 


r 

e3su'L 

1st  term  on 

2nd  line  of  (1-6) 

[(N-1)/N]g£. 

total  2nd 
line  of  (1-6) 

0.4 

.083 

-.147 

.030 

-0.117 

1.04 

.245 

-.250 

.276 

0.026 

1 

2.08 

-.173 

-.080 

-.025 

-0.105 

3.04 

-.389 

-.019 

-.034 

-0.053 

3.68 

-.426 

-.008 

-.017 

-0.025 

■ 

Table 

1-1  gives 

an  indication  of  the 

size  of  the 

"extra"  terms 

relative  to  Ejs 

u ^ for  some  selected 

r values. 

We  see  from  the 

table 

that  the 

extra  terms  are  small 

for  r > 3.0 

but  are  not  small 

in  general;  thus,  it  is  not  clear  that  the  BM  procedure  is  valid. 
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APPENDIX  J 


COMMENTS  ON  TRANSLATING  U vs  V TO  G vs  P 

The  important  point  is  that  getting  the  slope  (P  = -dU/dV)  from 

published  curves  of  U vs  V is  difficult  to  do  accurately.  People  who 

do  this  frequently,  feel  that,  even  with  the  typical  3"  x 4"  figure  from 

a journal  article,  it  is  better  to  assign  numerical  x,y  values  from  the 

figure  and  run  them  through  a spline  routine  than  to  try  to  take  the 

30 

slope  from  the  figure  directly. 


TABLE  J-l.  DETERMINATION  OF  P,G  FROM  STEVENSON'S  METALLIC 


U vs  V PLOT  (REF.  i) 


Volume 

i units)  (BM  units) 

pa 

(Mbar) 

PV 

(ev) 

U 

(ev) 

G 

(ev) 

13.15 

21.85 

1.27 

17.35 

0.11 

17.46 

13.93 

23.14 

0.984 

14.20 

-0.73 

13.47 

15.0 

24.92 

0.659 

10.25 

-1.63 

8.62 

15.8 

26.25 

0.494 

8.09 

-2.10 

5.99 

16.7 

27.70 

0.384 

6.66 

-2.50 

4.16 

17.5 

29.07 

0.321 

5.84 

-2.80 

3.04 

20.0 

33.22 

0.192 

3.98 

-3.44 

0.54 

25.0 

41.53 

0.106 

2.76 

-4.15 

-1.39 

30.0 

49.83 

0.080 

2.48 

-4.65 

-2.17 

35.0 

58.14 

0.0406 

1.48 

-4.96 

-3.48 

40.0 

66.44 

0.0218 

0.910 

-5.104 

-4.194 

50.0 

83.06 

0.0144 

0.745 

-5.300 

-4.555 

57.1 

94.85 

0.0 

0.0 

-5.360 

-5.360 

aFrom  DSPLSW  6/28/77 


*D.  J.  Stevenson,  Nature  258,  222  (1975). 

TQ 

Ray  Scanlon,  private  communication. 


12 

The  spline  routine  we  use  , DSPLSW,  fits  a cubic  function  over 
segments  (preserving  continuity).  Our  first  attempts  with  DSPLSW  led 
to  non-monotonically  decreasing  (magnitude)  slopes.  A better  choice  of 
end  points  and  end-point  slopes  seemed  to  improve  matters. 


TABLE  J-2.  DETERMINATION  OF  P,G  FROM  BM'S  METALLIC 


U vs  V PLOT  (REF.  2) 


Volume 
(BM  units) 

Pa 

(Mbar) 

PV 

(ev) 

U 

(ev) 

G 

(ev) 

13.8 

1.65 

14.20 

0.28 

14.48 

15.0 

1.35 

12.65 

-1.10 

11.55 

17.88 

0.76 

8.55 

-2.96 

5.59 

18.95 

0.60 

7.10 

-3.42 

3.68 

20.0 

0.47 

5.90 

-3.76 

2.14 

23.52 

0.24 

3.54 

-4.51 

-0.97 

26.85 

0.167 

2.69 

-4.92 

-2.23 

30.0 

0.121 

2.27 

-5.21 

-2.94 

35.0 

0.067 

1.45 

-5.48 

-4.03 

40.0 

0.048 

1.20 

-5.66 

-4.46 

49.0 

0.014 

0.430 

-5.846 

-5.42 

55.0 

-0.010° 

-0.327 

-5.850 

60.0 

-0.018? 

-0.684 

-5.804 

70.0 

-0.024° 

-1.029 

-5.670 

fFrom  DSPLSW  6/22/77 
“Computed  pressure  is  negative 


2M.  J.  M.  Bernal  and  H.  S.  W.  Massey,  Mon.  Not.  R.  Astr.  Soc.  114, 
172  (1954). 

l2We  are  indebted  to  Ray  Scanlon  for  furnishing  this  program  and  for 
discussions  related  to  its  use. 
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!j 

) 

, Selected  values  of  the  pertinent  entities  are  given  in  Tables  J-l 

i through  J-4  for  Stevenson  metallic,  BM  metallic,  piesent  metallic,  and 

Stevenson  (Ross)  mixture,  respectively.  For  the  volume  entries  the 
Stevenson  (S)  units  are  cm3/mol  NH4;  the  BM  units  are  (cm3/N  atom)  x 
1024.  The  U values  in  Tables  J-2  and  J-3  differ  somewhat  from  those  in 
columns  three  and  four  of  Table  5 due  to  smoothing  of  the  U vs  V curves 
by  the  DSPLSW  routine;  this  is  more  noticeable  at  the  low  volume  (high 
pressure)  end. 


TABLE  J-3.  DETERMINATION  OF  P,G  FROM  THE  PRESENT 


U vs  V CALCULATION 


Volume  Pa  PV  U G 

(BM  units) (Mbar) (ev) (ev) (ev) 


10.414 

1.905 

12.40 

1.50 

13.90 

14.83 

0.895 

8.29 

-2.07 

6.62 

17.45 

0.523 

5.67' 

-3.21 

2.46 

20.40 

0.253 

3.22 

-3.90 

- 0.68 

23.55 

0.132 

1.935 

-4.253 

- 2.32 

30.93 

0.074 

1.415 

-4.702 

- 3.29 

35.20 

0.057 

1.260 

-4.876 

- 3.62 

39.72 

0.041 

1.013 

-5.015 

- 4.00 

55.81 

0.010 

0.342 

-5.228 

- 4.89 

68.64 

0.004 

0.161 

-5.284 

- 5.12 

76.63 

-0.000 

0.000 

-5.293 

- 5.29 

83.31 

-0.001b 

-0.071 

-5.290 

100.0 

-0. 002b 

-0.120 

-5.271 

aFrom  DSPLSW 

U 

6/24/77 

'Computed  pressure  is  negative 
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TABLE  J-4.  DETERMINATION  OF  P,G  FROM  THE  STEVENSON  (ROSS)  MIXTURE 
U vs  V PLOT  (REF-  1) 


Volume  Pa  PV  U G 


(S  units) 

(BM  units) 

(Mbar) 

(ev) 

(ev) 

(ev) 

8.4 

14.0 

0.895 

7.81 

-4.01 

3.80 

8.94 

14.85 

0.802 

7.43 

-4.50 

2.93 

10.0 

16.6 

0.533 

5.55 

-5.24 

0.31 

11.25 

18.7 

0.374 

4.35 

-5.79 

-1.44 

12.5 

20.75 

0.319 

4.14 

-6.26 

-2.12 

15.0 

24.92 

0.124 

1.925 

-6.777 

-4.85 

17.5 

29.07 

0.0741 

1.345 

-7.026 

-5.68 

20.0 

33.22 

0.0527 

1.092 

-7.186 

-6.09 

25.0 

41.53 

0.0212 

0.550 

-7.384 

-6.83 

30.0 

49.83 

0.0068 

0.210 

-7.435 

-7.23 

aFrom  DSPLSW  6/27/77 

*D.  J.  Stevenson,  Nature  258,  222  (1975). 
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Units  Conversion: 

atomic  units  *-*■  cm3 

4-rr  3 r0. 52917A  10_8cm_  3 

Volume  = -lii  (a  in  a.u.)3  [ • — 6 ] ; 

3 a.u.  A 

3 

thus,  volume  (in  -c-m — ) = 0.62069  x 10"24  (a  in  a.u.)3, 

N-atom 


or,  volume  (in  x 1024)  = 0.62069  (a  in  a.u.)3. 


or,  volumeJJ.  = (a  in  a.u>)3. 

0.62069 

Stevenson  (S)  volume  units  +-*■  Bernal -Massey  (BM)  volume  units 


volume  = 10-24  (BM  volume  #)  — — — 

N-atom  ’ 


volume  = (S  volume  #) 


mol  NH,* 


10-24  (BM  vol  #)  — g-1"3  [ 6^02x1 02 3 N-atoms-j 

N-atom  mol  Nlh 


= 0.602  (BM  vol  #) 


mol  NH/i 


thus,  the  volume  in  S units  = 0.602  (volume  in  BM  units) 


dU/dV  <-+  P 

dU  = A(BM  U ev/N-atom) [-1.602xl0~12  erg-j 

dV  A(BM  V x 10- 2 4 cm3 /N-atom)  ev 

s.  rl  dyne-cm-i  r 1 Mbar  ■, 
erg  101 2dyne/cm2 

= (BM  slope  value) (1.602)  Mbar; 

thus,  (1.602) (slope  as  from  BM)  gives  P in  Mbar,  and  (0.602) (1 .602) 
'X (slope  as  from  Stevenson)  gives  P in  Mbar. 
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